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INTRODUCTION 


The  horseshoe-vortex  and  associated  corner  flows  which  occur  when  a 
blunt  obstruction  is  placed  within  an  approaching  boundary  layer  represent  a 
fundamental  three-dimensional  viscous  flow  of  considerable  interest  and 
importance.  Examples  of  this  type  of  flow  include  the  flows  near  an  aircraft 
wing/fuselage  junction  and  near  a  submarine  hull/sail  junction.  The  feature 
common  to  all  horseshoe  (or  necklace)  vortex  flows  is  that  a  non-uniform 
velocity  is  an  approaching  boundary  layer  meets  a  local  region  of  adverse 
pressure  gradient  due  to  the  blockage  effect  of  the  obstruction.  This  causes 
a  three-dimensional  boundary  layer  separation  and  the  formation  of  one  or 
more  horseshoe  vortices  around  the  obstruction.  In  addition  to  the  leading 
edge  region,  the  associated  corner  flows  downstream  of  the  leading  edge  are 
also  of  interest,  since  they  contain  streamwise  vortices  which  affect  the 
performance  of  both  the  airfoil  or  strut  and  also  other  devices  located 
downstream. 

The  problem  of  horseshoe-vortex/corner  flow  has  been  investigated 
previously  [1-2]  by  numerical  solution  of  the  Navier-Stokes  equations. 
Solutions  for  laminar  flow  past  an  elliptical  leading  edge  mounted  normal  to 
a  flat  plate  endwall  have  been  computed  for  both  zero  and  five-degree  angles 
of  incidence,  with  chordal  Reynolds  number  of  400,  and  Mach  number  of  0.2. 
Turbulent  flow  cases  have  been  computed  for  both  unswept  and  45-degree  swept 
elliptical  leading  edges  mounted  on  a  flat  plate,  with  Reynolds  number  of 
310,000  and  Mach  number  of  0.05.  Results  from  these  flow  calculations  have 
been  reported  by  Briley  and  McDonald  [1-2].  Calculations  for  a  blunt-fin 
induced  shock  wave  and  boundary-layer  interaction  flow  containing  a  horseshoe 
vortex  have  been  reported  by  Hung  and  Kordulla  [3].  Detailed  experimental 
measurements  for  incompressible  turbulent  horseshoe  vortex  flows  have  been 
obtained  recently  by  McMahon,  Hubbartt  and  Kubendran  [4]  and  by  Moore  and 
Forlini  [5]. 

A  major  impediment  to  the  study  of  this  and  other  three-dimensional 
flows  by  solution  of  the  Navier-Stokes  equations  has  been  the  high  cost 
(computer  run  time)  of  computing  accurate  solutions.  The  high  cost  is 
attributable  to  the  large  number  of  grid  points  inherent  in  three  dimensions 
and  to  a  loss  in  convergence  rate  associated  with  the  use  of  locally  refined 
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(nonuniform)  grids  which  are  necessary  to  define  the  multiple  length  scales 
present  in  high  Reynolds  number  viscous  flows.  Another  difficulty  is  that 
accuracy  can  be  degraded  by  the  use  of  artificial  dissipation  terms  which  are 
added  to  central  difference  approximations  for  convective  terms.  Attempts  to 
reduce  this  source  of  error  by  using  small  amounts  of  dissipation  have  led  to 
instability  in  three  dimensions  using  the  present  algorithm.  The  present 
investigation  was  undertaken  to  acquire  an  improved  understanding  of  the 
stability  and  accuracy  of  the  three-dimensional  ADI  scheme  for  a  scalar 
convection  model  problem,  which  would  lead  to  improvements  in  accuracy  and 
efficiency  of  the  solution  algorithm  for  the  Navier-Stokes  equations. 

In  the  present  report ,  the  methods  used  in  solving  the  three-dimensional 
Navier-Stokes  equations  for  horseshoe  vortex  flows  are  first  reviewed.  The 
questions  of  stability  and  error  due  to  artificial  dissipation  are  then 
examined  for  scalar  convection  in  three  dimensions.  Finally,  solutions  for  a 
laminar  horseshoe  vortex  flow  are  computed  using  an  improved  algorithm  with 
reduced  dissipation. 


GOVERNING  EQUATIONS 

The  three-dimensional  compressible  Navier-Stokes  equations  are  solved 
here  for  low  Mach  number  and  with  an  assumption  of  constant  stagnation 
enthalpy.  For  these  conditions,  steady  flow  solutions  closely  approximate  an 
incompressible  constant  density  flow  (cf.  Briley,  McDonald  and 
Shamroth  [6]).  A  zonal  approach  is  used  wherein  the  flow  is  computed  only  in 
a  subregion  of  the  overall  flow  field,  near  the  leading  edge  (Fig.  1).  The 
form  of  the  governing  equations  solved  permits  the  use  of  general 
nonorthogonal  body-fitted  coordinate  systems,  and  is  obtained  by  a 
transformation  from  Cartesian  to  general  nonorthogonal  coordinates.  The 
Cartesian  velocity  components  u{  and  density  p  are  retained  as  dependent 
variables  in  the  transformed  system  of  equations.  The  pressure,  temperature 
and  stagnation  enthalpy  are  denoted  p,  T  and  hQ,  respectively. 

All  variables  are  nondimens ional  in  the  present  formulation,  having  been 
normalized  by  reference  quantities  denoted  by  a  subscript  'r'.  The 
quantities  pr,  Ur ,  Tr  amd  Lr  denote  reference  values  for  density, 
velocity,  temperature  and  length,  respectively.  The  reference  pressure, 
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enthalpy  and  time  are  taken  as  PfUr  .  cpTr  and  Lr/Ur,  respectively, 

where  Cp  is  the  specific  heat  at  constant  pressure.  The  specific  heat 

ratio  is  Y,  and  Mr  =  Ur/cr  is  a  reference  Mach  number,  where  cr  is 

.  2 

the  reference  sound  speed  defined  by  cr  =  yRTr ,  and  R  is  the  gas 
constant.  The  reference  Reynolds  number  Re  is  defined  by  prUrLr/ur, 
where  ur  is  a  reference  viscosity. 

The  transformation  T  from  Cartesian  coordinates  x£  to  computational 
coordinates  yJ  is  given  by 


T  -  yj(xi)  i,j  -  1,2,3 

Spatial  derivatives  are  transformed  according  to 


(1) 


(2) 


where  unless  otherwise  stated  the  summation  convention  is  used  for  repeated 

indices,  and  y^ .  =  3yJ/3x..  The  coordinate  system  is  defined  by  specifying 

,11 

the  Cartesian  coordinates  of  each  computational  grid  point.  The  partial 
derivatives  3x./3y^  of  the  inverse  transformation  T  *  =  x-(yJ)  are  then 
computed  using  three-point  second-order  difference  formulas  with  uniform 
spacing  of  the  computational  coordinates  yJ  .  For  convenience,  the  yJ 
coordinates  are  normalized  to  give  a  unit  mesh  spacing  AyJ  =  1  for  each 
coordinate.  The  transformation  derivatives  3yJ/3x^  are  then  computed 
from  3x^/3yJ  using  standard  procedures  for  computing  derivatives  of 
inverse  functions  (cf.  Kaplan  (7]). 

The  transformed  Navier-Stokes  equations  can  be  written  in  the  following 
nondimensional  form:  The  continuity  equation  is 


3p 

3t 


+ 


0 


(3) 
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The  kth  component  of  the  momentum  equation  is  given  by 
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where  is  the  Kronecker  delta  function.  The  shear  stress  is  given 

by 
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The  equation  of  state  and  definition  of  stagnation  enthalpy  can  be  expressed 
for  a  perfect  gas  as 


p  -  PT/YM^  (6) 

hQ  =  T/(y-l)Mr2  +  q2/2  (7) 

where  q2  =  61Ju^uj.  Although  it  is  not  necessary,  it  is  both  convenient 
and  computationally  worthwhile  for  the  present  problem  to  assume  that  hQ  is 
a  constant  and  to  omit  solution  of  the  energy  equation.  This  results  in 
negligible  error  for  steady  flow  at  low  Mach  number  with  no  heat  addition. 
Equations  (6)  and  (7)  can  then  be  combined  to  produce  an  adiabatic  equation 
of  state 


p  =  P (h  -  q2/2)  (y-D/y 
o 


(8) 


which  is  used  to  eliminate  pressure  as  a  dependent  variable  in  Eq .  (4). 


METHOD  OF  SOLUTION 


The  basic  algorithm  considered  here  has  been  described  by  Briley  and 
McDonald  [8,  9]  and  employs  a  formal  time  linearization  to  produce  a 
noniterative  fully-coupled  approximation  for  nonlinear  systems  of  equations, 
which  is  solved  in  block-implicit  form  using  an  ADI  scheme  with  consistent 
intermediate  steps.  For  a  linear  scalar  diffusion  equation,  this  algorithm 
reduces  to  a  classical  ADI  scheme  considered  by  Douglas  and  Gunn  [10]. 

Warming  and  Beam  [11,  12]  have  introduced  a  very  concise  derivation  of  this 
same  algorithm  using  approximate  factorization  of  the  linearized 
approximation  written  in  'delta'  form.  The  works  of  Pulliam  and  Steger  [13], 
Thomas  and  Lombard  [14]  and  Shamroth,  McDonald  and  Briley  [15]  are 
representative  of  numerous  investigations  which  have  employed  this  basic 
algorithm. 


Linearization  and  Time  Differencing 

The  nonlinear  system  of  governing  equations  is  first  written  (at  a 
single  grid  point)  in  the  following  form: 


3H(<j>)  / 3t  =  D(<J>)  +  S(4>) 


where  <J>  is  the  column-vector  of  dependent  variables,  H  and  S  are 
column-vector  algebraic  functions  of  $,  and  D  is  a  column  vector  whose 
elements  are  the  spatial  differential  operators  which  generate  all  spatial 
derivatives  appearing  in  the  governing  equation  associated  with  that  element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit 
time-difference  approximation  of  (9): 

(Hn+1  -  Hn)/At  =  6(Dn+1  +  Sn+1)  +  (1-6)  (Dn  +  Sn)  HO) 

where,  for  example,  Hn+^  denotes  H(4>n+^)  and  At  =  tn+^  -  tn.  The 

parameter  8  (0.5  <  8  <_  l)  permits  a  variable  time-centering  of  the  scheme, 

2 

with  a  truncation  error  of  order  [At  ,  (8  -  1/2  At)]. 
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A  local  time  linearization  (Taylor  expansion  about  of  requisite 

formal  accuracy  is  introduced,  and  this  serves  to  define  a  linear 
differential  operator  L  such  that 

Dn+1  -  Dn  +  Ln  (4>n+1  -  4>n)  +  0  (At2)  (11a) 


Similarly, 


Hn+1  -  Hn  +  (3H/3$)n  (t1*1  -  *“)  +  0  (At2)  (lib) 

Sn+1  -  sn  +  (3S/3«>)n  (<j>n+1  -  4>n)  +  0  (At2)  (11c) 

Equations  (lla-c)  are  inserted  into  Eq .  (10)  to  obtain  the  following  system 
which  is  linear  in  <j>n+l 


(A  -  8At  Ln)  ($n+1  -  *n)  -  At  (Dn  +  Sn)  (12) 

and  which  is  termed  a  linearized  block  implicit  (LBl)  scheme.  Here,  A 
denotes  a  square  matrix  defined  by 

A  =  (3H/3$)n  -  6At  (3S/3$)n  (13) 

Equation  (12)  has  0  (At)  accuracy  unless  H  =  $,  in  which  case  the  accuracy  is 
the  same  as  Eq .  (10). 

Special  Treatment  of  Diffusive  Terms 

Spatial  cross-derivatives  are  present  in  viscous  terms  and  in  added 
artificial  dissipation  terms  of  the  present  formulation,  and  these  cross 
derivative  terms  are  evaluated  explicitly  at  tn.  To  preserve  notational 
simplicity,  it  is  understood  that  all  cross-derivative  terms  appearing  in 
Ln  are  neglected  but  are  retained  in  Dn .  In  addition,  although  diffusion 
coefficients  in  viscous  and  dissipation  terms  are  generally  functions  of  the 
dependent  variables,  these  coefficients  are  not  linearized  and  instead  are 
evaluated  implicitly  at  tn  during  each  time  step.  Notat ional ly ,  this  is 
equivalent  to  neglecting  derivatives  of  these  coefficients  with  respect  to 
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in  Ln,  which  are  formally  present  in  the  Taylor  expansion  (11a),  but 
otherwise  retaining  all  terras  in  both  Ln  and  Dn . 

It  is  important  to  note  that  neglecting  terms  in  Ln  has  no  effect  on 
steady  solutions  of  Eq.  (12),  since  —  4>n  =  0  and  thus  Eq.  (12) 

reduces  to  the  steady  form  of  the  equations:  Dn  +  Sn  =  0.  Aside  from 
stability  considerations,  the  only  effect  of  neglecting  terms  in  Ln  is  to 
introduce  an  0(At)  truncation  error. 

Consistent  Splitting  of  the  LSI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (12)  is  split 
using  ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional 
operator  L  is  rewritten  as  the  sum  of  three  "one-dimensional"  sub-operators 
L.£  (i  =  1,  2,  3)  each  of  which  contains  all  terms  having  derivatives  with 
respect  to  the  i-th  spatial  coordinate.  The  split  form  of  Eq .  (12)  can  be 
derived  either  by  following  the  procedure  described  by  Douglas  and  Gunn  in 
their  generalization  and  unification  of  scalar  ADI  schemes,  or  using 
approximate  factorization.  In  either  case,  for  the  present  system  of 
equations  the  split  algorithm  is  given  by 


(A  - 

BAtL*) 

(**  - 

$n)  =  At  (Dn  + 

sn) 

(14a) 

(A  - 

BAtL^) 

** 

<♦  ■ 

-  4>n)  =  A  ($*  - 

4>n) 

(14b) 

(A  - 

BAtL^) 

(*n+1 

-  <fn)  =  A  ($** 

-  4>n) 

( 14c ) 

where  <p*  and  <J>**  are  consistent  intermediate  solutions.  If  spatial 
derivatives  appearing  in  L,£  and  D  are  replaced  by  three-point  difference 
formuLas ,  then  each  step  in  Eqs.  (14a-c)  can  be  solved  by  a  block-tridiagonal 
e l iminat ion . 

Combining  Eqs.  (14a-c)  gives 


(A  -  8AtL")  A-1  (A  -  6AtL")  A_1  (A  -  6AtL^)  (4>n+1  -  $") 
=  At  (Dn  +  S°) 


(15) 
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which  approximates  the  unsplit  scheme  (12)  to  0  (At  ).  Since  the 
intermediate  steps  are  also  consistent  approximations  for  Eq.  (12),  physical 
boundary  conditions  can  be  used  for  and  $**.  Finally,  since  the  L{  are 
homogeneous  operators,  it  follows  from  Eqs .  (14a-c)  that  steady  solutions 
have  the  property  that  =  <p*  =  $**  =  $n  and  satisfy 

Dn  +  Sn  -  0  (16> 

The  steady  solution  thus  depends  only  on  the  spatial  difference  approximation 
used  for  (16),  and  does  not  depend  on  the  solution  algorithm  itself. 

STUDIES  ON  ACCURACY  AND  CONVERGENCE  RATE 

In  practical  high  Reynolds  number  calculations  using  the  split  LBI 
scheme  (14a-c),  artificial  dissipation  terms  are  added  to  the  central  spatial 
difference  approximations  for  convective  terms  to  control  numerical 
oscillations,  stabilize  the  solution  algorithm,  and  promote  convergence. 
Although  the  added  dissipation  terms  introduce  a  first-order  spatial 
truncation  error  in  steady  solutions,  the  magnitude  of  this  error  can  be 
controlled  using  adjustable  dissipation  parameters  and  is  at  worst  comparable 
to  that  associated  with  two-point  "upwind"  differences.  The  increase  in 
accuracy  derived  from  reduced  values  of  added  dissipation  has  been 
demonstrated  in  two-dimensions  by  Shamroth,  McDonald  and  Briley  [15]  in 
computations  of  viscous  transonic  flow  past  cascades  of  airfoils. 
Unfortunately,  instability  has  been  encountered  in  previous  three-dimensional 
flow  calculations  when  only  small  amounts  of  dissipation  are  used.  This  is 
not  surprising  since  the  Douglas-Gunn  ADI  (or  delta-form  approximate 
factored)  scheme  applied  to  a  simple  scalar  convection  equation  with  periodic 
boundary  conditions  is  known  to  lose  its  unconditional  stability  in  three 
dimensions.  This  three-dimensional  convective  instability  has  been  pointed 
out  by  Warming  and  Beam  [16],  Dwoyer  and  Thames  [17]  and  others.  In  the 
present  section,  the  behavior  of  the  ADI  scheme  (14a-c)  applied  to  scalar 
convection  in  three  dimensions  is  considered  further,  and  a  new  form  for  the 
artificial  dissipation  terms  is  examined. 


0 


0 


0 


0 


0 
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(constant  u)  is  approximated  on  a  unit  cube  with  an  equally  spaced  grid  such 
that 

*  JAx  ,  yk  *=  kAy  ,  zft  -  fcAz  ;  j.k.t  -  0,1,  ...N  (18) 

and  with  Ax  =  Ay  =  Az  =  h  =  N-*.  The  notation  Q1}  denotes 

J  >  K  ,  * 

<J>(x_jj  Yk»  zl >  tn)>  and  one  or  more  of  these  indices  will  often  be 
omitted  for  simplicity.  Equation  (17)  is  approximated  by  the  following  ADI 
scheme  analogous  to  (14a-c): 

(1  -  BAt  D  )  A4>*  «  At  (D  +  D  +  D  )  <j>n 
x  y  x  y  z 

**  * 

(1  -  BAt  Dy)  A4>  =  A<|> 

(1  -  BAt  D  )  A4>  =  A<J>** 
z 

where  A<f>  =  4>n+^  -  <J>n  and  At  =  tn+^  -  tn.  Central  spatial  differences 

are  used,  so  that  the  spatial  difference  operators  appearing  in  (19a-c)  are 

given  by 

Dx  =  (u/2h)  (<frj+1  -  ♦j_1)  (20 

with  analogous  expressions  for  Dy  and  Dz. 


(19a) 

(19b) 

(19c) 
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To  examine  the  stability  of  (19a-c)  using  the  von  Neumann  method,  a 
discrete  Fourier-component  solution  of  the  form 


rn,-,  i(wlxj  +  w2yk  +  W3Z£) 

k  £(w)  -  5  (w)  e  J 


is  substituted  into  (19a-c)  after  eliminating  the  intermediate  solutions  A<j>* 
and  and  the  amplification  factor  C  (At,  w)  =  ^  is  determined. 

Here,  u>  denotes  the  vector  of  frequencies  (o>i,  W2>  0)3),  and  i  =  -1.  The 
algorithm  is  stable  for  a  given  time  step  if  for  all  frequencies  w  present. 

+  S  “  Ul  <  1  (22) 

To  cover  all  possible  frequencies,  the  algorithm  must  be  stable  for  arbitrary 
10  in  the  range 

-  n/h  <_  -  n/h  (23) 


and  for  these  conditions,  it  has  been  found  that  the  algorithm  (19)  is 
unconditionally  unstable.  However,  this  particular  stability  analysis 
applies  rigorously  only  for  an  infinite  domain  with  periodic  initial 
conditions,  and  does  not  account  for  the  influence  of  nonperiodic  boundary 
conditions  on  a  finite  domain. 

In  the  present  investigation,  the  computation  of  steady  solutions  of  the 
Navier-Scokes  equations  is  of  primary  interest,  and  this  invariably  leads  to 
the  specification  of  nonperiodic  boundary  conditions  in  one  or  more 
directions.  In  addressing  steady  solutions  of  the  scalar  convection 
equation  (17),  it  is  natural  to  specify  function  values  at  inflow  boundaries 
and  to  employ  some  form  of  extrapolation  or  one-sided  differencing  at  outflow 
boundaries.  Accordingly,  the  behavior  of  the  algorithm  is  investigated  for 
u  >  0  and  with  the  (implicit)  boundary  conditions 
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{x  =  o 
y  *  0 
z«0 

.2  n+1  n  . 

6  A  -  0  at  x  -  1 
x 

,2, n+1  _  . 

6  $  *0  at  y  “  1 

y 


(24a) 

(24b) 

(24c) 


(24d) 


(25) 


2  2 

with  analogous  definitions  for  6^  and  6^.  The  extrapolation  outflow 

condition  is  equivalent  to  replacing  the  central  difference  formula  for  the 

direction  normal  to  the  outflow  boundary  by  two-point,  one-sided, 

differences  at  points  adjacent  to  the  outflow  boundary. 

The  stability  of  the  algorithm  (19)  subject  to  the  implicit  boundary 

conditions  (24)  was  tested  in  numerical  experiments  using  3  =  1.0  and 

different  choices  of  At,  h,  and  initial  conditions.  By  chosing  ^  =  0  on 

the  inflow  boundaries,  the  exact  solution  of  the  steady  difference  equations 

is  <j>  =0,  and  consequently,  the  value  of  ♦*!  ,  .  also  represents  the  error 

j  ,  k , 

e?  .  from  the  steady  solution.  The  stability  and  degree  of  convergence  was 
j  ,  k ,  Jo 

assessed  by  observing  the  behavior  of  en  with  increasing  n.  The  L2  norm  of 
the  error  at  interior  grid  points 


(26) 
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was  monitored  as  an  indicator  of  both  stability  and  rate  of  convergence  to 
the  steady  solution. 

In  each  of  the  cases  tested,  the  algorithm  (19)  with  boundary  conditions 
(24)  and  with  0  =  1.0  was  found  to  be  stable  for  sufficiently  small  Courant 
number  C  =  uAt/h,  in  contrast  to  the  unconditional  instability  indicated  by 
the  Fourier  analysis  with  frequencies  given  by  (23).  For  example,  Fig.  2 
shows  the  computed  error  behavior  for  50  steps  using  C  =  0.8  and  the 
(discontinuous)  initial  condition 


G(a)  = 


4>  -  G(x)  G(y)  G(z) 

2a  for  0  _<  a  <_  0.5 

l-2a  for  0.5  <  a  <  1 


(27a) 


(27b) 


The  results  in  Fig.  2  indicate  a  stable  calculation  with  good  error  reduction 
for  both  N  =  10  and  N  =  50.  In  another  example,  the  amplification  factor  for 
C  =  1.02  and  N  =  10  has  a  theoretical  value  of  1.00484  (unstable)  for  the 
frequency  u)  =  (-tt,  -2n,  5x).  Using  C  =  1.02,  N  =  10  and  the  initial 
condition 

<)>  »  G(-irx)  G(-2wy)  G(5irz)  (27c) 

G(a)  *  sin(a)  +  cos(a)  (27d) 

which  contains  only  this  single  frequency  and  does  not  satisfy  the  boundary 
conditions,  the  L2  error  was  reduced  to  less  than  10“^  its  initial  value  in 
100  iterations.  Had  the  theoreteical  amplification  factor  of  1.00484  been 
accurate  for  these  boundary  conditions,  this  error  would  have  increased  by 
62%  in  100  iterations.  The  discrepancy  between  the  instability  predicted  by 
the  (periodic)  Fourier  analysis  and  the  empirically  observed  stability  is 
attributed  to  the  non-periodic  boundary  conditions  (24).  Since  these 
boundary  conditions  are  treated  implicitly,  they  influence  the  solution 
simultaneously  at  all  grid  points  during  each  time  step. 

A  heuristic  treatment  of  non-periodic  boundary  conditions  within  the 
von  Neuman  approach  can  be  accomplished  by  taking  the  region 
(0  <_  x,  y,  z  <_  1)  to  be  a  half  interval  for  each  coordinate,  with  periodic 
extension  of  the  (nonperiodic)  solution  into  the  full  interval 
(-1  <_  x,  y,  z,  <_  1).  At  interior  grid  points  (1  £  j,  k,  l  <_  N-l),  the  error 
ej,k,i.  from  the  steady  solution  can  be  expressed  as  a  finite  Fourier  series 
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containing  only  sine  terms  (odd  extension),  with  frequencies  given  by 


us. 

-  m,  it 

m- 

»  1,2,. 

.  .N-l 

1 

1 

1 

(28a) 

“2 

-  m2ir 

m2 

-  1,2,. 

.  .N-l 

(28b) 

-  m.ir 

-  1,2,. 

.  .N-l 

(28c) 

This  representation  correctly  satisfies  the  condition  of  zero  error  at  inflow 
boundaries,  but  incorrectly  assumes  zero  error  at  outflow  boundaries.  The 
amplification  factor  C(At,  us)  for  the  restricted  range  of  frequencies  (28) 
provides  a  useful  estimate  of  the  observed  stability,  although  it  is 
heuristic  in  its  approach  to  the  non-periodic  boundary  conditions.  A 
numerical  computation  for  the  frequency  range  (28)  with  N  =  10  indicates  that 
the  maximum  |c|  numbers  for  which  ojj  =  a>2  =  103  =  us,  and  therefore  it  is 
sufficient  to  consider  the  behavior  of  C(At,  us)  for 


0  <  <uh  <  ir 


(29) 


For  6  =  J,  this  leads  to  the  following  stability  condition  (See  Appendix  A) 

C  -  uAt/h  <_  /3/2  (30) 

The  stability  of  the  algorithm  (19)  with  boundary  conditions  (24)  was 
tested  in  numerical  experiments  for  N  =  10  with  initial  conditions  given  by 


♦  K  S 

*  t-  c= 


sin(jTtx)  sin(kiry)  sin(Hrtz) 


(31) 


which  includes  contributions  of  unit  amplitude  from  all  frequencies  in  (28). 
The  quantity 


is  taken  as  an  empirical  measure  of  the  maximum  amplification  factor  \{i  for 
the  purpose  of  assessing  stability.  Computed  values  of  <j)  after  30  iterations 
are  shown  in  Fig.  3,  along  with  theoretical  curves  for  both  the  full  range  of 
frequencies  (23)  and  for  positive  frequencies  in  the  range 

0  _<  (d^,  (^2 ,  —  "/h  (33) 

which  led  to  the  estimate  (30)  for  stability.  The  theoretical  curves  were 
obtained  numerically  using  a  sampling  increment  of  n/lOh  for  each  frequency. 
For  the  conditions  of  these  test  calculations,  the  stability  condition  (30) 
provides  a  useful  and  somewhat  conservative  estimate  of  the  observed  limit  of 
stability.  Although  the  present  results  provide  only  a  limited  model  for  the 
Navier-Stokes  equations,  the  stability  and  rapid  convergence  rates  obtained 
here  with  proper  time  step  selection  are  very  encouraging.  It  is  later 
demonstrated  here  that  rapid  convergence  can  also  be  obtained  for  the 
three-dimensional  Navier-Stokes  equations  with  small  values  of  artificial 
dissipation. 

Second-Order  Artificial  Dissipation 

One  method  of  adding  dissipation  is  equivalent  to  a  replacement  of  the 
convective  derivative  operator  for  each  coordinate  direction  by  a  modified 
operator,  as  follows: 

3  .  3  I  u  1  Ax  3  ( 34 ) 

3x  3x  °  2  .2 

3x 

with  analogous  replacements  for  v3(  )/3y  and  w3(  )/3z.  When  three-point 
central  differences  are  used  on  a  uniform  grid,  the  resulting  approximation 
is  equivalent  to  the  two-point  upwind  difference  scheme  when  0=1.  This 
latter  scheme  has  first-order  accuracy  and  is  especially  inaccurate  when  the 
velocity  is  not  aligned  with  the  computational  grid.  Alternatives  to  this 
upwind  scheme  have  received  considerable  attention  (see,  for  example, 

Raithby  [18],  Baliga  and  Patankar  [19],  and  Brooks  and  Hughes  [20]).  A  very 
simple  but  effective  method  of  improving  the  accuracy  of  the  upwind  scheme  is 
to  express  the  scheme  in  terms  of  central  differences  and  an  added 
dissipation  as  in  (34)  and  then  reduce  the  dissipation  parameter  o.  For 
sufficiently  small  values  of  o,  this  method  obviously  approaches  the  accuracy 


of  the  second-order  central  difference  scheme.  Dissipation  of  the  form  given 
in  (34)  will  be  referred  to  here  as  'upwind'  dissipation. 

When  the  velocity  is  not  aligned  with  the  coordinate  mesh,  the  accuracy 
for  a  given  value  of  a  can  be  improved  by  adopting  a  different  form  for  the 
dissipation  terms,  which  acts  to  smooth  convected  quantities  only  in  the 
direction  along  streamlines.  This  form  of  dissiaption  is  used  in  the 
streamline-upwind  method  of  Brooks  and  Hughes  [20]  and  is  also  present  in  the 
tensor  viscosity  method  of  Dukowicz  and  Ramshaw  [15].  Dissipation  of  this 
type  will  be  referred  to  here  as  'streamwise'  dissipation,  and  can  be 
illustrated  as  follows: 

If  s  denotes  distance  along  a  streamline,  and  if  s  is  the  unit  vector  in 
the  direction  of  the  streamline,  then  the  velocity  vector  U  can  be  expressed 
as  U  =  qs  where  q2  =  U  •  U,  and  the  convective  operator  can  be  expressed  as 

U-V-qf^  (35) 


where 


3s 


s  •  7 


(36) 


Streamwise  dissipation  can  be  added  by  replacing  the  convective  operator  by  a 
modified  operator  analogous  to  (34)  as  follows: 


U  •  V  t- 


qAs  3 


(37) 


where 


(qAs)^  =  (u  Ax)^  +  (v  Ay)^  +  (w  Az)^ 


(38) 


In  three  dimensions  and  for  constant  velocity,  the  streamwise  dissipation 


The  above  form  is  equivalent  to  upwind  dissipation  only  in  one  dimension 
(u  *  0,  v  =  w  =  0)  or  if  the  mixed-derivative  terms  are  omitted.  In  the 
numerical  algorithm,  the  mixed  derivative  terras  are  evaluated  explicitly. 
Warming  and  Beam  [22]  have  shown  for  a  three-dimensional  scalar  diffusion 
equation  that  the  explicit  treatment  of  mixed-derivative  terms  does  not  upset 
the  unconditional  stability  of  the  ADI  scheme  (19)  with  0=1. 

The  effect  of  upwind  and  streamwise  dissipation  on  accuracy  was  examined 
for  a  three-dimensional  test  problem  suggested  by  Abarbanel,  Dwoyer  and 
Gottlieb  [23].  Here,  the  convection  equation  (17)  is  solved  in  the  unit  cube 
with  equal  mesh  increments  and  for  N  =  10.  An  exact  steady  solution  $  of 
equation  (17)  for  these  conditions  is  given  by 

4>  =  cos  [2n(x-y)]  +  cos  [2ir(y-z>]  +  cos  [2ir(x-z)]  (40) 

In  this  test  problem,  the  velocity  vector  is  directed  diagonally  to  each 

incremental  mesh  cube.  In  the  present  calculations,  the  correct  boundary 

conditions  for  <j>  were  prescribed  from  this  known  exact  solution  ( <j>  =  $  at 
2  2 

inflow;  <5  <J>  =  6  ♦  at  outflow).  It  was  noted  that  the  central  difference 
scheme  without  dissipation  (o  =  0)  reproduced  the  exact  solution  at 
all  grid  points.  Consequently,  the  error  <f>  -  $  in  the  present  calculations 
is  entirely  due  to  the  added  dissipation.  The  maximum  error 
|<{>  -  ^/(^max  -  $min)  at  interior  points  is  shown  in  Fig.  4  as  a  function 
of  the  dissipation  parameter  a  for  both  the  upwind  (34)  and  streamwise  (39) 
forms  of  dissipation.  Although  the  coarse  mesh  used  (N  =  10)  causes 
relatively  large  errors  at  large  values  of  o,  the  streamwise  dissipation  (39) 
provides  a  considerable  improvement  in  accuracy  for  a  given  value  of  o. 

Other  solutions  were  computed  in  which  resolution  was  improved  by  using  the 
same  coarse  grid  (N  =  10)  but  for  computational  domains  smaller  than  the  unit 
cube.  These  results  confirm  that  the  increased  accuracy  of  streamwise 
dissipation  over  upwind  dissipation  for  fixed  o  is  maintained  when  the  grid 
resolution  is  increased. 
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COMPUTED  RESULTS  FOR  A  LAMINAR  HORSESHOE  VORTEX  FLOW 


Flow  Conditions 


Solutions  are  presented  here  for  laminar  flow  at  zero  incidence  past  an 
elliptical  leading  edge  geometry  mounted  between  parallel  flat  plate 
endwalls.  The  purpose  of  these  calculations  was  to  determine  whether  the 
improved  understanding  of  stability  and  accuracy  acquired  in  the  present 
model  problem  studies  could  be  exploited  to  allow  computation  of  accurate 
solutions  of  the  three-dimensional  Navier-Stokes  equations  using  the  split 
LBI  scheme  with  small  amounts  of  artificial  dissipation.  The  leading-edge 
geometry  for  the  present  calculations  is  the  same  as  that  considered  by 
McMahon,  Hubbartt  and  Kubendran  [4]  for  turbulent  flow  conditions.  A  laminar 
flow  case  was  chosen  here  to  allow  an  accurate  assessment  of  the  numerical 
method,  which  is  not  clouded  by  extraneous  factors  associated  with  turbulence 
modelling.  The  computation  of  turbulent  flow  cases  will  be  undertaken  in  a 
future  investigation. 

The  flow  geometry  (Fig.  1)  consists  of  a  strut  of  constant  thickness  W 
having  an  elliptical  leading  edge  with  1.5:1  ratio  of  major  to  minor  axis. 

The  strut  is  mounted  normal  to  parallel  flat  plate  endwalls  whose  separation 
distance  is  5.0W,  and  whose  leading  edges  are  located  a  distance  6.0W 
upstream  of  the  leading  edge  of  the  strut.  The  length  L  of  the  strut  within 
the  computational  domain  is  2.5W.  The  flow  considered  has  a  Reynolds  number 
Re  =  200  (based  on  L)  and  Mach  number  Mr  =  0.1,  each  based  on  upstream  flow 
conditions . 

Boundary  Conditions 

Since  the  computational  domain  is  chosen  to  be  a  region  in  the  immediate 
vicinity  of  the  leading-edge/corner  flow  geometry  (cf.  Fig.  1)  embedded 
within  a  larger  overall  flow  system,  inflow  and  outflow  boundary  conditions 
which  adequately  model  the  interface  between  the  computed  flow  and  the 
remainder  of  the  flow  system  are  required.  The  inf  low/ out  flow  conditions 
used  are  derived  from  an  assumed  flow  structure  and  are  chosen  to  provide 
inflow  with  prescribed  stagnation  pressure  (and  stagnation  enthalpy)  in  an 
inviscid  core  region  and  with  a  given  axial  velocity  profile  shape  in  the 


endwall  boundary  layer,  and  to  provide  outflow  with  a  prescribed  distribution 
of  static  pressure  in  the  cross  section.  This  approach  to  inflow/outf low 
boundary  conditions  has  been  discussed  previously  in  [1-2]  and  is  only 
summarized  here. 

First,  the  boundary  layer  thickness  <5(xi)  on  the  endwall  flat  plate  is 
approximated  by  its  distribution  from  the  Blasius  flat  plate  solution.  At 
the  inflow  boundary,  a  "two-layer"  boundary  condition  is  employed  such  that 
stagnation  pressure  pD  is  fixed  at  the  free  stream  reference  value  in  the 
core  flow  region  (y3  >  5)  and  an  axial  velocity  profile  shape 

3  3 

ul/ue  =  f (y  /S)  is  fixed  within  the  boundary  layer  region  (y  <_  &) .  Here, 
ue  is  the  local  edge  velocity  which  varies  with  time  and  is  adjusted  after 
each  time  step  to  the  value  consistent  with  pQ  and  the  local  edge  static 
pressure,  which  is  determined  as  part  of  the  solution.  The  remaining  inflow 
conditions  are  U2  =  3  03/ 3n  =  3  p/ 3n  =0,  where  n  denotes  the  normal 
computational  coordinate,  y1.  For  outflow  conditions,  a  constant  static 
pressure  is  imposed,  and  second  derivatives  of  each  velocity  component  are 
set  to  zero.  At  no-slip  surfaces,  each  velocity  component  u{  is  set  to 
zero,  and  the  remaining  condition  applied  to  these  surfaces  is  that  the 
derivative  of  pressure  in  the  direction  normal  to  the  surface  is  zero.  This 
condition  approximates  the  normal  momentum  equation  to  order  Re-1  for  viscous 
flow  at  a  no-slip  surface.  The  final  boundary  to  be  considered  is  the  plane 
parallel  to  the  endwall  and  in  the  free  stream.  This  boundary  is  assumed  to 
be  a  plane  of  symmetry,  so  that  the  flow  represented  is  that  past  the  strut 
mounted  between  parallel  flat  plates . 

Computational  Details 


For  constant  stagnation  enthalpy,  the  governing  equations  consist  of 
continuity  (3)  and  three  components  of  momentum  (4).  After  eliminating 
pressure  using  (8),  these  equations  are  solved  using  the  split  LBI 
scheme  (14),  linearized  with  respect  to  the  dependent  variables 
<(>T  =  (p,  ui  ,  U2 ,  U3).  Central  spatial  differences  are  used  on  an  equally 
spaced  grid  for  each  transformed  coordinate  yJ ,  and  streamwise  dissipation 
terms  analogous  to  (39)  are  added  to  these  equations. 

For  each  governing  equation,  the  streamwise  dissipation  term  can  be 
expressed  as 
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(41) 


2<l£  si  8* 

2  ayV 


where  sJ  =  UJ/Q  are  components  of  a  unit  vector,  and  where 

uJ  “  y]±  ui  (42a) 

Q2  =  UjUk  6jk  (42b) 

Here,  o  is  a  dissipation  parameter,  ^  is  a  scalar  representing  the 
appropriate  dependent  variable  (p  for  continuity  and  u^  for  the  kth 
momentum  equation),  and  k  is  0  for  continuity  but  otherwise  1. 


The  quantity  V  defined  by 


V 


where 


i  2 
^  (y\) 

i  ’x 


(43) 


(44) 


is  a  measure  of  the  physical  diffusion  in  the  streamwise  direction.  The 
artificial  dissipation  term  (41)  is  omitted  at  points  for  which  V  is  greater 
than  op^Q/2.  In  the  present  calculations,  the  value  o  =  0.1  was  used  for 
the  dissipation  parameter.  It  should  be  recalled  that  for  scalar  convection 
in  one  dimension,  the  value  o  =  1.0  is  equivalent  to  two-point  upwind 
differencing,  and  thus  the  value  o  =  0.1  is  relatively  small. 

Although  a  uniform  grid  was  used  for  the  transformed  coordinates  yJ , 
the  grid  used  is  highly  nonuniform  in  physical  coordinates  x£,  as  shown  in 
Fig.  I,  and  was  chosen  to  provide  resolution  of  several  length  scales  known 
to  be  present  for  this  type  of  flow.  Care  was  taken  to  provide  resolution  of 
the  boundary  layers  on  the  strut  and  endwall  surfaces,  of  the  shear  layer 
near  the  leading-edge  stagnation  line,  and  of  the  corner  flow  region  very 
near  the  endwal 1/ leading-edge  intersection.  Two  computat ional  grids, 
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consisting  of  15  and  29  equalLy  spaced  points  for  each  transformed  coordinate 
direction,  were  used  to  compute  solutions  which  were  otherwise  identical. 

For  the  finer  grid,  the  mesh  spacings  adjacent  to  the  endwall  and  strut 
surfaces  were  0.0005W  and  0.0085W,  respectively. 

A  pseudo-time  dependent  iteration  procedure  was  employed  in  which  the 
square  matrix  A  in  the  LBI  scheme  (14)  is  replaced  by  a  modified  matrix 
A  *  AB,  where  B  is  a  diagonal  conditioning  matrix.  This  technique  was 
employed  by  Briley,  McDonald  and  Shamroth  [6]  to  improve  the  convergence  rate 
of  this  algorithm  when  applied  to  low  Mach  number  flows.  In  the  present 
calculations,  a  diagonal  matrix  B  whose  diagonal  elements  are 
(l/YMr  ,  1,  1,  1)  was  used.  In  addition,  a  spatially  varying  time  step  was 
used  to  avoid  instability  analogous  to  that  encountered  for  the 
three-dimensional  scalar  convection  equation,  and  to  improve  the  convergence 
rate.  Several  procedures  for  time  step  selection  for  the  Navier-Stokes 
equations  are  presently  under  investigation,  and  the  present  results  should 
be  taken  ony  as  a  demonstration  that  this  technique  can  provide  a  good  rate 
of  convergence  to  the  steady  solution.  A  more  detailed  account  of  variuos 
techniques  for  time  step  selection  will  be  given  in  a  forthcoming  paper. 

Computed  Results 

Before  proceeding  to  a  discussion  of  the  horseshoe  vortex  flow,  an 
indication  is  given  of  the  degree  and  rate  of  convergence  obtained  for  these 
three-dimensional  calculations,  and  of  the  effect  of  halving  the  mesh  on  the 
computed  solutions.  The  convergence  rate  for  each  of  the  two  grids  is  shown 
in  Fig.  5,  where  the  maximum  residual  in  the  field  is  normalized  by  its  value 
at  the  start  of  the  calculation.  A  very  good  rate  of  convergence  was 
obtained  for  a  relatively  small  value  of  dissipation  parameter  (o  =  0.1). 
Reducing  the  normalized  residuals  to  10“  gave  convergence  to  within  a 
tolerance  not  distinguishable  in  the  plots  of  computed  flow  variables  given 
here,  and  this  required  about  30  to  60  iterations,  respectively,  for  the 
(15  x  15  x  15)  and  (29  x  29  x29)  grids. 

The  effect  of  halving  the  mesh  spacing  on  the  computed  distributions  of 
flow  variables  at  a  representative  location  is  shown  in  Fig.  6.  The 
distributions  of  each  velocity  component  and  static  pressure  coefficient 
Cp  =  2(p  -  l/YMr2)  in  the  X3  direction  are  given  for  a  line  normal  to  the 
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endwalL  and  within  the  plane  of  symmetry  upstream  of  the  strut  leading  edge, 
located  a  distance  0.22W  upstream  of  the  leading  edge.  This  location  has  the 
peak  value  of  reversed  flow  in  the  uj  velocity  component.  The  comparison  in 
Fig.  6  is  representative  of  that  found  throughout  the  flow  field.  The 
solution  does  not  change  significantly  when  the  mesh  is  halved,  and  this  is 
an  indication  that  the  solutions  are  accurate  and  that  the  flow  behavior  is 
adequately  resolved  using  these  grids.  The  sensitivity  of  the  coarse  mesh 
solution  to  the  location  of  the  outflow  boundary  was  also  tested  by  adding 
grid  points  to  move  the  outflow  boundary  downstream  a  distance  6.1W,  and  this 
had  no  significant  affect  on  the  solution. 

Computed  results  for  the  (29  x  29  x  29)  grid  are  shown  in  a  series  of 
plots  in  Figs.  7,  8  and  9a-f.  In  Fig.  7,  velocity  vector  plots  in  the  plane 
of  grid  points  adjacent  to  the  no-slip  endwall  surface  are  contrasted  with 
corresponding  results  from  the  symmetry  plane  midway  between  the  endwalls. 

The  vector  magnitudes  are  renormalized  for  each  plot  and  thus  indicate  only 
flow  direction  and  relative  magnitude  within  the  plot.  A  reversed  flow 
region  within  the  horseshoe  vortex  is  clearly  visible  near  the  endwall.  The 
reversed  flow  region  includes  two  saddle-point  flow  separations  associated 
with  the  horseshoe  vortex  flow  structure.  One  of  the  saddle  points  is  very 
close  to  the  leading  edge  and  is  not  clearly  visible  in  Fig.  7. 

The  computed  flow  behavior  in  the  corner  region  near  the  intersection  of 
the  strut  and  endwall  is  shown  in  Figs.  8-9.  The  flow  is  shown  in  planes 
normal  to  the  strut  surface  and  intersecting  the  endwall,  as  indicated  in 
Fig.  8.  Vector  plots  of  velocity  within  these  planes  and  contour  plots  of 
normal  velocity,  static  pressure  coefficient  Cp  and  total  pressure 
coefficient  CpQ  are  shown  in  Figs.  9a-f.  The  cross  sections  have  a  width 
of  2W  and  height  of  1.25W.  The  horseshoe  vortex  formation  in  the  stagnation 
symmetry  plane  (plane  l)  can  be  seen  in  Fig.  9a.  There  is  a  strong  downward 
velocity  toward  the  endwall  in  the  region  near  the  leading  edge  (behind  the 
saddle-point  separation),  with  maximum  downward  velocity  about  45  per  cent  of 
the  free  stream  reference  velocity.  The  maximum  reversed  flow  velocity  (uj) 
is  about  7  per  cent  of  the  reference  velocity  (see  also  Fig.  6).  The 
horseshoe  vortex  continues  its  development  at  downstream  locations  and  moves 
outward  away  from  the  corner  region  and  into  the  free  stream.  The  contours 
of  velocity  normal  to  the  cross  sectional  planes  give  an  indication  of  the 
shear  layer  development  on  the  strut  and  endwall  surfaces  and  also  show  the 
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distortion  of  the  shear  layers  in  the  corner  region  due  to  the  horseshoe 
vortex.  The  static  pressure  field  is  highly  three-dimensional,  with 
significant  spanwise  gradients  along  the  strut  surface  and  low  pressure  near 
the  strut/endwall  intersection.  Although  no  other  analytical  results  or 
experimental  measurements  of  the  three-dimensional  horseshoe  vortex  flow  are 
available  for  comparison,  the  present  computed  results  are  consistent  with 
flow  visualization  studies  of  related  leading  edge  vortex  flows. 

CONCLUDING  REMARKS 

Significant  progress  has  been  made  in  improving  the 
accuracy/convergence-rate  properties  of  the  present  numerical  method  for  the 
three-dimensional  Navier-Stokes  equations.  It  has  been  demonstrated  for  the 
laminar  flow  considered  that  an  accurate  (grid  insensitive)  steady  solution 
can  be  computed  and  that  good  convergence  rates  can  be  obtained  in  three 
dimensions,  using  only  small  amounts  of  artificial  dissipation.  These 
calcultions  were  performed  using  locally-refined  nonuniform  grids  of  the  type 
needed  to  define  the  multiple  length  scales  present  in  high  Reynolds  number 
viscous  flows.  A  future  study  will  address  the  computation  of  turbulent  flow 
cases  with  resolution  of  the  viscous  sublayer  region.  Based  on  our  recent 
experience  with  a  vectorized  version  of  this  code  on  a  CRAY-1  computer,  it 
should  soon  be  possible  to  compute  solutions  for  a  30  x  30  x  30  grid  (27,000 
points)  in  about  10  minutes  of  run  time. 


APPENDIX  A 


Under  the  condition  that  the  frequencies  for  each  coordinate  direction 
are  equal  (wi  =0) 2  =  w 3  =  <u) ,  the  amplification  factor  for  the  algorithm  given 
in  (19a-c)  can  be  written  in  the  form 


“l  +  1  h 

r  m  - - — — 

(  A.  1 ) 

a2  +  1  62 

where 

-  a2  ■  1  -  3  (CA) 2 

(A. 2a) 

82  -  -  3CA  +  (CA)3 

(A. 2b) 

8l  =  82  +  3CA/6 

(A. 2c) 

and  where 

C  *  uAt/h 

(A. 3) 

A  *  Bsin(a)h) 

(A. 4) 

Because  ai  =  02,  neutral 

stability  occurs  when  81  = 

02 »  which  yields 

CA  =  [3(26-l)/26]1/2 

(A. 5) 

The  maximum  value  of  A  occurs  for  oah  -  ^/2,  and  this 

gives  the  following 

stability  condition: 

C  <  8-1  [3(26-1) / 2 S ] 1 ^ 2 

(A.  6) 

The  special  case  8  =  1  gives  C  *3/2 
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Fig.  9c  -  Computed  Flow  Variable.';  in  Plane  3. 
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Fig.  9d  -  Computed  Flow  Variables  in  Plane  4. 
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Fig.  9e  -  Computed  Flow  Variables  in  Plane  5. 
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Fig.  9f  -  Computed  Flow  Variables  in  Plane  6. 
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